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Abstract. 

Properties of reaction zones resulting from A + B — > C type reaction-diffusion 
processes are investigated by analytical and numerical methods. The reagents A and B 
are separated initially and, in addition, there is an initial macroscopic inhomogeneity 
in the distribution of the B species. For simple two-dimensional geometries, exact 
analytical results are presented for the time-evolution of the geometric shape of the 
front. We also show using cellular automata simulations that the fluctuations can be 
neglected both in the shape and in the width of the front. 
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1. Introduction 

Front propagation has attracted considerable research activity during the past years. 
A large number of phenomena in different fields (e.g. crystal growth, chemical pattern 
formation, and biological invasion problems) are determined by the properties of various 
types of fronts. 

An important class of problems is related to the case of front propagation into 
unstable states |JJ. They describe situations in which the dynamics of a system is 
governed by a front which invades a domain where the system is in an unstable state, 
leading to the onset of a new stable state (often with spatio-temporal structures present). 
Examples include population dynamics modeled by the extended Fisher-Kolmogorov 
equation [2], onset of vortices in Taylor-Couette flows |3] and of convection rolls in 
Rayleigh-Benard experiments (modeled by the Swift-Hohenberg equation |lj), phase 
turbulence in autocatalytic chemical reactions modeled by the Kuramoto-Sivashinsky 
equation [Hj. 

Most of the studies have been addressing the properties of reaction-diffusion 
fronts in homogeneous media [Hj. Presently, however, an increasing importance 
is attached to the more complicated cases of fronts propagating in inhomogeneous 
media. Inhomogeneities may have different natures: chemical (related to spatial 
inhomogeneities in the concentration of the reagents), physical (reagents' diffusion 
coefficients may vary in space), geometrical (changes in the geometry of the medium), 
etc. Also, inhomogeneities may be either related to an intrinsic disorder present in 
the system, or may be created artificially. There have been both theoretical and 
experimental effort in studying various situations, see e.g. [7J |HJ HE EH IUJ • Examples 
are numerous in the context of biology [12], population dynamics jTHj, heterogeneous 
catalysis [JJ], or heterogeneous excitable media like the photosensitive Belousov- 
Zhabotinsky reaction ^Hl UHl Ej- In particular, it was argued that under some 
circumstances react ion- diffusion fronts may behave like light-waves when undergoing 
refraction and reflection [9] ITU] ITTj. 

In the present paper we shall be concerned with a different class of phenomena 
which emerges when the unstable state is produced by reaction processes in a front 
moving diffusively. In the wake of the front, the dynamics of the unstable state may 
generate various spatio-temporal patterns. A prototype of this pattern-forming process 
is the Liesegang phenomena [IE]t where a chemical reactant (called inner electrolyte) 
is dissolved in a gel matrix and a second reactant (outer electrolyte) diffuses in and 
reacts with the inner electrolyte. Under certain conditions, the dynamics of the reaction 
product generates a family of high-density precipitate zones whose properties obey 
generic laws [20] which have been understood in terms of a simple phase separation 
scenario [21] • 

Most of the Liesegang pattern studies have been made in homogeneous 

| While the reagents in Liesegang phenomena are chemicals, examples where the reagents are of more 
exotic nature (such as topological defects ^Hj) are a l so known. 



Reaction- diffusion fronts with inhomogeneous initial conditions 



3 



conditions [TBI, when the initial concentration of the inner electrolyte is uniform. In 
recent work by Grzybowski et al [221, however, the experiments were performed in 
inhomogeneous conditions (with a step-like profile in the height of the gel containing 
the inner electrolyte). It was argued that Liesegang bands emerging in the two regions 
of different thicknesses join each other at the step according to a Snell law, where the 
indices of refraction are related to the so-called spacing coefficients of the patterns. 

In order to develop a theoretical basis for the above findings, we shall 
first investigate the properties of the diffusive reaction front in the presence of 
inhomogeneities in the initial distribution of the inner electrolyte. This is the 
simplest way to take into account the inhomogeneities and it allows to obtain detailed 
understanding and analytical solutions. The ultimate goal, left for further studies, 
will be the use of the reaction zone results as an input to the description of the band 
formation in the spinodal decomposition scenario [21] which has been shown to correctly 
describe the pattern formation in the homogeneous case. 

The paper in organized as follows. In Sec. 2, we present a brief review of the 
main results known about the properties of diffusive reaction fronts in case when 
the electrolytes are initially separated but otherwise are homogeneously distributed. 
Section 3 is devoted to the study of fronts when inhomogeneity is present in the initial 
state of the inner electrolyte. Both analytical results, at the mean-field level, as well as 
numerical results are given. Cellular automata simulations are presented in Sec. 4 and 
the role of fluctuations is discussed. 

2. Reaction front: Homogeneous distribution of the inner electrolyte 

We review here the properties of reaction front emerging in an A + B — > C process 
where the particles A and B diffuse and react with rate k upon contact. Initially, the 
reagents are separated, with the inner electrolyte B occupying the x > half-space 
(homogeneously distributed with concentration bo), while the outer electrolyte A is 
homogeneously distributed with concentration a in the x < half-space. Usually, the 
case ao bo is considered and one finds that A and B react in a confined region (the 
reaction front) which moves in the positive x direction. The properties of this front 
can be addressed at different levels: a mean-field like one, neglecting fluctuations, or 
approaches that take the fluctuations into account. 

2.1. Mean- field analysis 

The mean-field study of the front resulting from the reaction A + B — > C has been made 
by Galfi and Racz [23]. The generalization to the reaction nA + mB —>■ C has been 
given by Cornell et al |24) . We shall recall here the main results without going into the 
details of the calculations. 

In the homogeneous initial condition case and with the fluctuations neglected, the 
concentration of the reagents a(x, t) and b(x, t) will depend only on the variable x at all 
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times and the problem becomes one-dimensional. The reaction-diffusion equations for 
the concentrations a(x, t) and b(x, t) read: 

d t a = D a d 2 a — kab , d t b = D b d 2 b — kab , (1) 

where k is the reaction rate and D a , D b are, respectively, the diffusion coefficients of A 
and B. 

The production rate of C, which is a quantity of interest since it is the source 
for the precipitation process in Liesegang phenomena |21j . is obtained by computing 
R(x,t) = ka(x , t)b(x , t) . Actually, the region where R(x,t) is significantly different 
from zero defines the reaction zone. The center of the zone, Xf(t), can be found from 
the condition a(xf,t) = b(xf, t), and for the case of D a = Db = D this leads to especially 
simple results. Indeed, one notices then that u(x, t) = a(x,t) — b(x,t) obeys a simple 
diffusion equation. Introducing q = b /a , the solution is given in terms of the error 
function erf, 



u(x, t) — a , , 

2 2 \2^/Dt / 

and the condition a(xf,t) = b(xf,t) yields 



1 — q 1 + q „ / x 
erf 



(2) 



This means that the front moves diffusively and its diffusion coefficient Df is given by 



the equation eri(^jD f /2D) = (1 - q)/(l + q). 

In order to complete the description of the front one needs the solution of the 
nonlinear part of the reaction-diffusion equations, 

d t a = Dd^a — k(a 2 + au) . (4) 

It is not possible to obtain a general solution of this equation. However, assuming that 
the width w(t) of the reaction front, defined as 

/oo 
(x — Xf) 2 R(x, t)dx 

= -^r* , (5) 

R(x, t)dx 



is much smaller that the diffusion length ^iffW = v^Dt (an assumption that is justified 
a posteriori), one finds that, in the long-time limit, the quantities of interest take the 
following scaling forms in the reaction zone region: 

w(t)~t a , a{x,t) = t~^a{xr a ) , 

b(x, t) = r^b(xr a ) , R(x, t) = t-PR{xr a ) . (6) 

Here the exponents are given as a = 1/6, (3 = 2/3 and 7 = 1/3. The exponents actually 
obey two general scaling relations following from the facts that ^diff(t) w(t) and that 
the reaction zone is fed by the fluxes of diffusing particles [23 121] , 

a + 7 = l/2, 2a + 7 = /3. (7) 
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In the mean-field approximation, one has an extra scaling relation, namely (3 = 27, 
and this leads to the above values of the exponents. Note that the case D a ^ is 
more complicated. However, the exponents obey the same scaling relations, and only 
the scaling functions a, b, and R are affected. 

A similar analysis can be carried out for the more general reaction nA + mB — ► 
C [21] ■ The main modification arises when expressing the production rate, and one has 
now the scaling relation (3 = (n + m)j, leading to the mean-field exponents: 
n + m — 1 n + m 1 

a = T( — , — rr\> p = — , — rr> 7 = ^ — rr- ( 8 ) 

2(n + m + 1) n + m + 1 n + m + 1 

I?. I?. Beyond mean- field 

Numerical simulations have been performed for the general reaction nA + mB — > C [24J, 
for several values of (n, m), including the case n — m = 1, and confirm the validity of the 
above scaling relations. Moreover, the mean-field exponents are recovered in dimensions 
larger than d = 2. In d — 1, the exponent values strongly differ from their mean-field 
values. This leads to the hypothesis that the upper critical dimension for diffusive 
reaction fronts is d u — 2. 

The validity of this hypothesis has been supported by theoretical arguments 
(Cornell and Droz |25J) based on the scaling behavior of the steady-state front generated 
by imposing antiparallel fluxes of A and B particles. This situation is much easier to 
investigate, since the front is time-independent, and depends only on three relevant 
parameters, the flux J, the diffusion coefficient D, and the reaction rate k. The results, 
however, can be directly applied to the time-dependent case too, since the front is formed 
quasistatically by diminishing fluxes J it) ~ t~ x l 2 . Assuming that scaling applies even 
when fluctuations are present (as seen in numerical simulations), a dimensional analysis 
and scaling arguments lead to the following conclusions: first, there is an upper critical 
dimension d u = 2/(m + n — 1) above which mean-field predictions are correct; second, 
for d < d u (only realizable for m = n = 1), the width exponent is a = l/2(d+ 1), which 
for d — 1 gives a = 1/4 instead of a = 1/6 for the mean-field case. 

In summary, one finds that, in the physically relevant dimensions, the reaction 
zone moves diffusively and its width is negligible with respect to the diffusion length 
of the problem. There is a third result which has not been discussed here but will be 
important in a later application when the front will figure as a source of particles in 
the description of the formation of Liesegang bands. Namely, the front leaves behind a 
constant concentration of C particles [2lH l2Hj. 

3. Reaction front: Inhomogeneous distribution of the inner electrolyte 

There are many ways to introduce inhomogeneities. We shall consider here a simple 
case which, on one hand, can be treated analytically and, on the other hand, mimics a 
system with a geometrical inhomogeneity quite similar to the experimental setup used 
in [22] ■ Namely, we shall assume that the inhomogeneity is in the initial concentration 
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of the inner electrolyte. As described in Fig^ the reagents are separated in the plane 
(x, y) into three zones, E , Si, and E 2 . For x < and Vy, a = a ; b = (region E ), 
while for x > the half-plane is divided in two domains, separated by a line making an 
angle 9 with the x axis. In the upper domain Ex, a = 0; b = &oi> while in the lower 
domain E2, a = 0; & = 602- This corresponds to an engineered initial inhomogeneity of 
the reagent B. 



a o 







X 



Figure 1. Geometry of the inhomogeneous media considered. 

As before we restrict ourselves to the reaction A + B — ► C, although the results 
could be extended to the more general case mA + nB — > C . We also assume that the 
diffusion coefficients are equal D a = D b = D. 

3.1. Mean- field description of the reaction front 

The mean-field description is again in terms of the reaction-diffusion equations (0), the 
only difference being that the concentrations depend on two spatial variables (x, y) and 
the diffusion terms are now expressed through the two-dimensional Laplacian. 

The position of the reaction front (xf(t) ,yf(t)) is again obtained by solving the 
diffusion equation for u(x, y, t) = a(x, y, t) — b(x, y, t) and finding (xf, yj) satisfying the 
condition u(xf, yf, t) = 0. The diffusion equation being linear, the solution u(x, t) for the 
initial condition in Fig|T]can be constructed as a superposition of two particular solutions 
u = U\ + u 2 . Here ui(x,y,t) satisfies the initial condition U\{x < 0,y,t = 0) = a and 
Ui(x > 0, y,t = 0) = —6 i and, according to Eq. (J2J), is given by 



Ui(x,y,t) = a 



1-qt 1 + qi , x 
— — erf 



2VDt 



(9) 



where q\ = b i/a . 
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The second function u 2 (x,y,t) satisfies the initial condition u 2 = in S and Si, 
and u 2 = &01 — b 2 = Qo(?i — 92) i n ^ 2 with q 2 = b 02 /a . The solution is constructed in 
terms of the Green function for the diffusion operator 

= (-^) < 10 > 

as an integral restricted to the region S 2 

«,(*,»,«) = *(*-*)/ *WG(«- (11) 

For simple cases, the integration can be performed analytically, leading to an 
explicit expression for u 2 . For example, for the symmetric case when 9 = 0, one obtains 

u(x,y,t) = 

{ * * \ZV L>1 J 

(12) 

The position of the reaction front is thus given by the implicit equation: 

1^ - X A^~ erf (X f ) + [1 + erf (X f )] [1 - erf (Y f )} = , (13) 

where 

X f = ^=, Y f =y$L. (14) 
J 2^Dt 1 2VDi 

are the scaled, time-independent coordinates of the front. The corresponding shape of 

the front is shown in the left panel of Fig El for different values of the parameters q\ and 

Far from the x axis (\y\ ^> y/2Dt), the reaction front is not affected by the 
inhomogeneity, it is parallel to the y axis and moves diffusively, along the x direction. 
The diffusion coefficients are given, respectively, by erf(yi5^ 2 /2D) = (1 — q 1 i2 )/(l+q , i, 2 ). 

Near the x axis, there is a crossover region between the two almost linear segments 
of the front. The crossover is smooth and can be characterized by the front position 
x°f(t) = 2Xf\^Dt on the x axis and by the angle ip between the tangent to the front at 
X° f and the x axis. The front position on the x axis is obtained from 

erf (X?) = 1 ~ fa + 92 |( 2 , (15) 
1 f) l + (?i + <&)/2' 1 ' 

and, of course, it also moves diffusively. The angle ip on the other hand is a constant 



[l + (gi + g 2 )/2] 



2 



tan ^ = — ^ ^" 1 exp[-(X» f y] . (16) 
Qi — 12 

A similar calculation can be performed for an arbritrary angle 9, however the results 
can no longer be expressed in terms of explicit functions. Namely, the position of the 
front (in terms of scaled variables) is given implicitly by the following equation: 
2-q 1 -q 2 2 + qi + q 2 



4 

Qi - <?2 r X 



er 



f(X f ) 



I f dzexp(-z 2 ) erf(Y) - X f t&n9 - zt&n9) = . (17) 

J — 00 



Reaction- diffusion fronts with inhomogeneous initial conditions 



8 




Figure 2. Shape of the reaction-diffusion front (in scaled time-independent 
variables) for different values of the parameters. Left panel: front shape for 
9 = 0. The thick continuous line corresponds to q± = 0.2 and qi = 0.02; the 
thin continuous line to q\ = 0.2 and q2 = 0.1; thick dashed line: q\ = 0.01, 
q2 = 0.001; and thin dashed-line: q\ = 0.01, qi = 0.005. X® gives the 
intersection with the x axis, and ip is the angle between the tangent at X9 
and the x axis. Right panel: front shape for different values of 9 (q± = 0.2 and 
q2 = 0.02). The thin straight lines indicate, respectively, the border between 
the regions Si and £2 for the considered values of 9. 

The profile of the front is given in the right panel of Fig. 2 for different values of the angle 
9. The analytical results have been checked by numerical integration of the reaction- 
diffusion equations for the reagents A and B. Note that the infinite reaction rate k limit 
is especially convenient for carrying out the numerical work since the position of the front 
does not depend on k for the D a = case. However, it does not allow the investigation 
of the width of the front. For finite k, the width of the front in the asymptotic regions 
Yf — > ±00 follows from the case of homogeneous initial distribution of inner electrolyte, 
i.e. w(t) ~ t a with a = 1/6. Investigations of the crossover region in the case when 
fluctuations are also included (see next section) indicate that no anomalous behavior 
occurs. Thus we expect that the width remains much smaller than the diffusion length 
and the scaling w(t) ~ t a with a = 1/6 remains valid everywhere along the front. 

The conclusions one can draw from the analytic results and from FigJ21 are as 
follows. First, the front motion remains diffusive but the shape deviates from straight 
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line. Second, the distortion is smooth and can be parameterized by the two asymptotic 
limits at Yf — > ±00, and by a linear parameterization of the crossover region. 

3.2. Beyond mean- field 

To test whether the mean-field predictions are robust against including fluctuations, we 
will compare the above results to those given by a mesoscopic approach, namely by the 
simulation of the reaction-diffusion front by a cellular automata. This approach also 
provides informations about how the width of the front scales with time, since a finite 
reaction rate is considered. 

The model (described in more detail in [26J) is defined on a two-dimensional square 
lattice of 250 x 250 sites with synchronous dynamics at discrete time steps. One iteration 
of the automata is defined as the combination of one diffusion and one reaction step. 
The diffusion coefficients for A and B particles are chosen to be equal, and the reaction 
A + B — > C takes place only for particles entering a site from opposite directions. We 
use thus a finite reaction rate, since the probability for a reaction of two particles is 
equal to 1/4. The front position is defined as the region where the difference of the two 
concentrations, computed in a Moore neighborhood for the A and B particles, vanishes. 
To allow for better statistics the results are averaged over 750 runs. 

The initial state is prepared as illustrated in FigO In the region £0 (with a length 
of 100 sites), all the sites are fully occupied, i.e. the concentration of A's is Oq = 1. The 
regions £1 and £2 (of length of 150 sites) are randomly occupied by B particles with 
concentration 6 01 = 0.2 and b Q2 = 0.02, respectively. 

The results of the cellular automata simulations for the case of 9 = are shown in 
FigOl We find good agreement between the analytical results and the averaged position 
of the front obtained from the cellular automata approach. The fluctuations in the 
system therefore appear to play no role for the position of the front. Equally good 
agreement was found for the case with a nonzero angle 9 = 20°. 

Preliminary investigations of the time dependence of the width of the front at the 
separation line between £1 and £2 seems to indicate the same time scaling w(t) ~ t 1 / 6 
as in the homogeneous case (as mentioned before, far from the inhomogeneity area the 
behavior is obviously the one following from homogeneous initial distribution of the 
inner electrolyte). 

4. Discussion 

We have considered the motion of the reaction zone for a case of simple initial-state 
inhomogeneity which mimics the experimental setup where the Snell law had been 
discussed for Liesegang bands [22]- Our main finding is that although the straight- 
line shape of the front becomes distorted due to the inhomogeneity, the basic properties 
such as (1) the diffusive overall motion, (2) the negligible width, and (3) the irrelevance 
of fluctuations in the physical dimensions remain unchanged. 
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250 




Figure 3. Results of the cellular automata simulations of the reaction-diffusion 
front for 9 = 0. The values of the initial concentrations are ao = 1, &oi = 0.2 
and &oi = 0.02, and the diffusion coefficients for the A and B particles are 
equal D a = Df,. The straight horizontal line separates the high- and the low- 
concentration regions of B particles. The points in the dotted region represent 
the instantaneous position of the front, calculated from condition u = a — b = 
after 3000 timesteps. The thin noisy line gives the position of the front averaged 
over 750 realizations of the process. The analytic result is shown by the thick 
line. 



The next step towards understanding whether a version of the Snell law is valid for 
Liesegang bands meeting at a line of inhomogeneity should be the use of our present 
results as an input for a theory of Liesegang band formation. The best candidate should 
be the spinodal decomposition description which has been shown [2T| to reproduce all 
the known generic properties of the Liesegang patterns. 
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